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Abstract 

The effects of annealing and irradiation on the evolution of Cu clusters in a-Fe are 
investigated using object kinetic Monte Carlo simulations. In our model, vacancies act as carriers 
for chemical species via thermally activated diffusion jumps, thus playing an important role in 
solute diffusion. At the end of the Cu cluster evolution, the simulations of the average radius and 
number density of the clusters are consistent with the experimental data, which indicates that the 
proposed simulation model is applicable and effective. For the simulation of the annealing process, 
it is found that the evolution of the cluster size roughly follows the 1/2 time power law with the 
increase in radius during the growth phase and the 1/3 time power law during the coarsening phase. 


In addition, the main difference between neutron and ion irradiation is the growth and evolution 


process of the copper-vacancy clusters. The aggregation of vacancy clusters under ion irradiation 
suppresses the migration and coarsening of the clusters, which ultimately leads to a smaller average 
radius of the copper clusters. Our proposed simulation model can supplement experimental 
analyses and provide a detailed evolution mechanism of vacancy-enhanced precipitation, thereby 
providing a foundation for other elemental precipitation research. 

1. Introduction 

The formation of solute-decorated clusters is one of the three main damage mechanisms in 
reactor pressure vessel (RPV) steels under irradiation, and this process is sensitive to the structure 
and composition of RPV steels. Generally, the alloying elements in RPV steels of A508-II mainly 
comprise C (0.19 wt.%), Mn (1.20-1.43 wt.%), Ni (0.73—-0.79 wt.%), Cr (0.06—0.12 wt.%), Mo 
(0.48—0.51 wt.%), and Cu (0.034—0.070 wt.%). In the past, the formation of Cu-rich clusters in 
iron has been widely studied experimentally [1-2]. Many studies have found that the formation of 
these clusters under service temperatures and neutron irradiation is a major cause of mechanical 
property degradation in RPV steels, such as hardening and embrittlement, which limits the lifespan 
of nuclear power plants. In earlier studies, it was observed that the precipitation of Cu solute in 
RPV steel could induce enrichment of the Mn/Ni/Si content; thus, the Cu content in RPV steel is 
strictly controlled. However, the influence of Cu on the evolution of matrix defects is of great 
significance for studies of the radiation-induced segregation process under irradiation; more 
importantly, it provides a reference for the subsequent study of the solute cluster precipitation of 
other elements. Therefore, vacancy-enhanced (Cu solute segregation affected by 
temperature/irradiation in Fe-Cu systems has often been studied. Using the positron annihilation 
technique, the enhanced diffusion of Cu was evidenced; meanwhile, vacancy dissociation from the 
Cu-vacancy cluster and the combination of clusters was observed [3-4]. Furthermore, Jin et al. 
investigated the effects of the irradiation dose and irradiation depth on the Cu cluster formation 
process using positron annihilation spectroscopy [5]. 

However, considering the difficulty of neutron irradiation experiments, computer simulations 
have been used as an alternative to neutron irradiation experimental studies [6-7], and these 
simulations can demonstrate the evolution of nanostructures under irradiation conditions and 
evaluate the mechanical properties of materials. Among common simulation methods, the Monte 
Carlo method is widely used for the analysis of reactor physics owing to its advantages in solving 


complex problems. Moreover, with the advancement of computer technology, its application for 


simulations of various physical reactor processes is increasing [8-10]. In fact, atomic kinetic Monte 
Carlo (AKMC) and object kinetic Monte Carlo (OKMC) methods allow large time-scale 
simulations over large spatial regions, up to the life span of a reactor vessel (30—40 y), and the 
microstructural evolution under irradiation in ferritic steels can be analyzed. 

AKMC is a type of rigid lattice model that can describe atomic motion, but crystallographic 
changes cannot be simulated. However, for OKMC, the atomic structures of the lattice and clusters 
are simplified and replaced by abstract objects, which are used to describe the properties of 
different types of defect structures through corresponding parameters. The AKMC model is based 
on the first neighbor jumps in a-Fe, thereby reproducing in detail the diffusion of defects or solutes, 
which leads to the formation of voids and precipitates. An introduction to the main features and 
different applications of the AKMC model is available in reference [11]. 

Correspondingly, OKMC simulations treat all defects as objects, including point defects (i.e., 
vacancies and self-interstitial atoms (SIAs)) and solute-decorated matrix defects (i.e., Cu-vacancy 
clusters). The characteristics of these objects, such as their shape, type, size, and position, are fully 
abstracted and parameterized. During the simulation, these objects follow several possible events, 
mainly migrating and dissociating or reacting with other objects, i.e., merging into larger clusters 
or annihilating. 

Soisson and Fu obtained AKMC simulation results that were consistent with experiments [12]. 
Considering the strong interaction between vacancies and copper clusters, their model simulated a 
series of vacancy transitions on the surface of copper clusters, which resulted in movement of the 
copper clusters. It was confirmed that the movement of clusters was an important mechanism for 
the precipitation of copper clusters in a-Fe. Furthermore, Castin and Pascute combined AKMC 
and OKMC simulations and also obtained good results consistent with those obtained 
experimentally: when the size of one copper cluster was larger than a certain threshold, the atomic 
structure was ignored, and the cluster simulation performed as an object [13]. This led them to 
extend the simulation scale to the coarsening state, allowing the simulation results to be compared 
with broader experimental results. 

Based on these previous studies, an OKMC simulation is proposed to describe the evolution 
of copper cluster segregation under irradiation in Fe-Cu systems. In this model, we further 


parameterize the properties of small clusters and apply a new method to define the factor for 


rescaling the simulation annealing time when dealing with the annealing process. This method 
enables us to simulate Cu segregation in a-Fe under irradiation. 

The simulation model in this study is built based on the framework of an open-source code 
(MMonCa) [14], and the fundamentals of our algorithm are described in the first part of this paper. 
Then, the evolution of Cu clusters in the Fe-Cu system during the annealing process is successfully 
simulated using the proposed model. Finally, we conduct a comparative study of the precipitation 
of Cu-vacancy clusters under neutron irradiation and ion irradiation, and the main evolutionary 
characteristics of these clusters are analyzed under different conditions. 

2. Object kinetic Monte Carlo simulation method 
2.1 OKMC model 

The kinetic Monte Carlo algorithm assumes that the system transitions randomly between 
different states. The transition is a Markovian process in which the transition probability, rj, 
depends only on the properties of the initial state, i, and the final state, j. The transition probability, 
rij, is the basic input parameter of the Monte Carlo algorithm. In this study, the transition probability, 
rj, is calculated using the Arrhenius formula [15]: 

rij = Pij X exp(— Eij/kpT), (1) 
where the pre-factor, Py, is the frequency of attempts, and Ey is the activation energy when the 
system transitions from state i to state j. After the transition probabilities, 7, of all possible states 
in the system are defined, we can perform an OKMC simulation using the residence time algorithm, 
which is the most commonly used algorithm in OKMC simulations [16,17]. 

2.1 Basic assumptions 

Generally, we assume that copper clusters migrate in the iron matrix through a vacancy 
mechanism, which implies a series of vacancy jumps at the surface of the copper clusters. However, 
because of the foundational assumptions of the OKMC simulation, the atomic structure of the 
cluster is simplified and the clusters are treated as independent objects with different properties, 
such as their size, trapping radius, attempt frequency, activation energy, and binding energy 
functions. Thus, a series of macro-equivalent descriptions of the properties of the cluster are 
needed, among which the two most important parameters are the diffusion coefficient and the 
migration coefficient; both of these parameters can be calculated from the AKMC simulation 
results of previous studies. For example, for the migration capacity of copper clusters, a vacancy 


and a pure copper cluster were placed in the simulation at a certain temperature, and the migration 


process of the vacancy in the lattice site was simulated; the calculation was stopped when the 
cluster produced an overall movement through the vacancy transition after a long period of 
operation. Using this process, a sufficient amount of data was obtained through repeated 
independent simulations, and the macroscopic equivalent diffusion coefficient of copper clusters 
could be obtained through statistical calculations [13]. The other basic parameters required in this 
study were also derived from the results of first-principles calculations and AKMC simulations in 
previous studies [13, 18, 19]. 

In addition, although the internal structures of defects and clusters are ignored, defect clusters 
in body-centered cubic (BCC) iron have two possible shapes in the OKMC simulation: spheres 
(3D) or dislocation loops (2D), where vacancy-type clusters usually exist as 3D spheres and SIA- 
type clusters usually exist as 2D dislocation loops. All objects move discretely in a given space 
according to their physical laws. For example, the mobile direction of the dislocation loops is along 
the direction of the Burgers vector. Moreover, defect clusters in the system interact during the 
evolution process to produce a defect reaction network with a corresponding reaction distance. 
When the distance between the two instances is less than the sum of their respective reaction radii, 
a reaction can occur. The reaction radius for each particle object and cluster object is defined as 
0.55 nm in this model. 

In particular, although there is a strong interaction between vacancies and Cu clusters, a stable 
Cu-SIA coupling is difficult to form in the Fe matrix; thus, copper atoms are considered to only 
interact with single vacancies or vacancy clusters but not with SIA clusters (dislocation loops) in 
this study. However, SIA clusters can still react with the vacancies in Cu-vacancy composite 
clusters, which would lead to elimination of the same number of vacancies in the Cu-vacancy 


composite cluster. The interactions of the defect clusters in this simulation are shown in Fig. 1. 


Dislocation loo 
Cu-vacancy pair e j 
“Gr ae Ete 
è SIA SOS oe? 
Vacanc 6 pe 
Y 8o? Ko) 
oe 3 
Cu-vacancy & ese 
cluster 
ove 


Fig. 1 (Color online) Reaction network of defect clusters in the Fe-Cu system 


Although the copper lattice has a face-centered cubic (FCC) structure, experiments have 
confirmed that when the diameter of the clusters is less than 5 nm, the copper clusters will maintain 
the same BCC structure as the iron matrix during their growth. When the diameter of the clusters 
is greater than 12 nm, the crystal structure of the copper clusters is transformed into the FCC 
structure [20, 21]. This conclusion has also been verified by molecular dynamics simulations, 
which further revealed that the existence of vacancy clusters in the Cu-vacancy composite cluster 
enhanced the structural stability of the cluster and inhibited the occurrence of this transition [21]. 
Overall, although the OKMC simplification ignores the description of FCC-structured copper 
clusters formed in BCC-structured iron, it is still acceptable to simulate the precipitation of BCC- 
structured copper clusters in iron, especially when the copper clusters are less than 5 nm in 
diameter. 

3. Parameterization 

In this study, the main parameters required for a cluster object include its mobility and stability, 
which are functions of the number of defects contained in the cluster. The abbreviations of these 
parameters are listed in Table 1, referring to vacancy-type clusters and SIA-type clusters, and the 
selection of their numerical values is described below. 


Table 1. Abbreviations of cluster object parameters 


Ns Number of defects contained in the cluster 
vê ig(Ns ) Attempt frequency for migration of the cluster 
Ep (Ns) Potential energy for defects to form in the cluster 
E¢(Ns) Formation energy of the cluster 
Ebina(Ns) Binding energy of the defects to the mother cluster 
E® (Ns ) Dissociation energy of a defect 


3.1 Vacancy-type clusters 

For a single vacancy, the migration energy has been calculated to be 0.64 eV [22] and 0.67 eV 
[23] by density functional theory (DFT); it has been experimentally determined to be 0.55-0.57 
eV [24], while the calculated value from the classical Mendelev potential is 0.63 eV [25]. In our 
study, the migration energy of the vacancy is set to 0.63 eV, between the above theoretical 
calculation and experimental values, to ensure good reliability. At the same time, we use the Debye 


frequency as the attempt frequency of the migration, v,;= 6x10!? s'. For vacancy clusters 


containing N, vacancies, when N, = 2—10, the AKMC simulation results [26] are used as the 
migration energies and corresponding transition attempt frequencies. For vacancy clusters with 
N, > 10, the calculated migration energy fluctuates between 1.15 and 1.50 eV, and thus an 
intermediate value of 1.30 eV is selected in this study. Similarly, it is assumed that the diffusion 
mechanism is surface diffusion [27,28]; that is, the movement of these vacancy clusters is the result 
of a series of surface vacancy transition events. Thus, the corresponding migration attempt 


frequency can be calculated as follows: 


V4 


vihig) = Gore (4) 

Before parameterizing the vacancy cluster stability, we must define the formation energy of 

vacancy clusters. In MMonCa, the cluster formation energy is the sum of the formation energies 
of the constituent particles minus the potential energy E,(N,) of the clusters, En(Ny) 

E¢(Ny) = Ny x Ef — Ep (Ny), (5) 
where Ef is the vacancy formation energy in the iron matrix, which is set as 2.02 eV in this study 
according to the value calculated in the AKMC simulation. Then, E,(N,) is the accumulation of 
binding energies for each unit vacancy bound to this cluster: 

Ep Ny) = Ep(Ny — 1) — Ebina (Ny — 1), (6) 
where Epina (Ny — 1) is the binding energy of a single vacancy to the cluster with N,-1. In this 
study, the corresponding AKMC results are used for the binding energies of vacancy clusters with 


N,= 2-4. For vacancy clusters with N, > 4, the fitting function in the literature [19] is used: 


EY,.a(N,) = 1.71 + 3.39716[(N,)? — (Ny + 13]. (7) 

The stability of vacancy clusters can then be defined to determine the rate of decomposition 

of vacancy clusters. Because the event in which a vacancy cluster emits another vacancy cluster 

has a very low probability [22], the main means of decomposing vacancy clusters is the continuous 

emission of single vacancies into the matrix. The dissociation energy, Eji,(Ny ), of vacancy 

clusters with N, is the sum of the binding energy of a single vacancy to the mother 
cluster, Eina (Ny ), and the single vacancy migration energy, Emig: 

Egis(Ny ) = Epinay ) + Enig- (8) 

It should be noted that the vacancy migration energy during dissociation from the mother 

cluster is equal to its migration energy when transitioning in the matrix, 0.63 eV, which is 


independent of the N, of the vacancy cluster. In addition, the dissociation attempt frequency of the 


cluster should be proportional to the size, N,, which is also the condition required for the cluster 
absorption and emission point defects to reach a balance. Then, the attempt frequency is defined 


in this study as follows: 


Viis (Ny ) = v1 X (Ny). (9) 

As mentioned previously, the movable copper clusters will contain at least one vacancy, 

whereas pure copper clusters without vacancies are regarded as immovable objects. Then, the Cu- 

vacancy clusters with N, = 1 are defined first. For Cu-vacancy clusters with Ncy= 2-5 and N, = 1, 

the frequency of transition attempts is based on the AKMC results [18]. For clusters with Ncy > 5 

and N, = 1, the mobility decreases with increasing Ncy following the surface diffusion mechanism 
[29]: 


vı 
(Ncu+Ny)4/3 ` o 


Vmig(Ncu + Ny) = 

At the same time, the migration energy of the Cu-vacancy pair is close to that of a single 
vacancy [19] according to DFT calculations, and thus it is also assigned a value of 0.63 eV. For 
clusters with Nc, = 2—4 and N, = 1, we use the AKMC simulation results mentioned previously 
[18]. For clusters with Ncy = 15-6000 and N, = 1, the migration energy is assigned a value of 1.2 
eV, as calculated in reference [13]. For clusters with Nc, = 5—14 and N, = 1, the corresponding 


value is taken as 1 eV. Figure 2 shows the decrease in Dy with increasing Ncy, which is similar to 


the AKMC results obtained in a previous study. 
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Fig. 2 Diffusion coefficient of Cu clusters at 500 °C 


For the stability of Cu-vacancy clusters with Ny = 1, only two types of dissociation events are 
considered in the simulation: the emission of single vacancies or Cu-vacancy pairs from their 
mother clusters. The corresponding dissociation energy is equal to the migration energy of the 
emitted defect object plus its binding energy to the mother cluster. For example, the dissociation 
energy, E};,, for a vacancy emission can be calculated as follows: 

Edis = Ebina (Ncw Nv ) + Emig» (11) 
where Exina(Ncw Ny ) is the binding energy of one vacancy to a cluster with Ncy and Ny; Emig is 
the migration energy of the vacancy, and the attempt frequency for emission is defined as 610!”. 


The migration energy of the emitted object is also considered to be independent of the size of its 


mother cluster, and Eyig= Pate vacancy = 0,63 eV for all Nc,,. We then calculate EYina (Necu Ny ) = 
Edis — Emig from the dissociation energies, Epj,. For clusters with Ncy = 2-5, the dissociation 
energies are obtained from reference [18]. For Ncy = 6, the values are extrapolated by fitting a 


function to the AKMC results using the following equation: 


Ena (Ncw 1) = 0.9 + 1.6 x (Neu)? — (Neu + 194). (12) 

This function is assumed to have the same functional form as Eq. (7), and the corresponding 
coefficients are then trained using the AKMC values. 

For Cu-vacancy clusters with Ny > 2, the cluster mobility also decreases as the number of 
vacancies increases, following the same surface diffusion mechanism described by Eq. (10). To 
determine the stability of these clusters, we use the results of molecular dynamics simulations that 
provide the binding energy of copper—vacancy composite clusters through computational 
simulations of metropolis Monte Carlo (MMC) and molecular dynamics binding [19, 29]; the 
stability is a binary function of N, and Ncy. In this study, the partial differential form is used as the 
input parameter, Eina (Necu Ny ). 

3.2 SIA-type clusters 

The method for defining the mobility of SIA-type clusters is similar to that for vacancy-type 
clusters, and the diffusion pre-factor is determined using the same equations. For SIA clusters of 
size N; < 7, the migration energy and migration attempt frequency used in this study are obtained 
from molecular dynamics studies [19, 29] and experimental studies [25] based on Mendeleev’s 
potential. It should also be noted that when the SIA cluster grows beyond a certain threshold, it 


will transform into 1/2<111> loops, and N; = 7 is selected as the threshold for this transition. At 


the same time, we do not distinguish between <111> dislocation loops and <100> dislocation loops 
when simulating their migration, and thus we use an effective migration energy and set a 
dislocation loop distinction threshold, Nj,. For dislocation loops with N; = 7 - Nip, it is assumed 
that only part are <100> dislocation loops. Therefore, we use the effective migration energy for 
the dislocation loops whose size is in this interval. It is assumed that the <111> dislocation loops 
are transformed into <100> dislocation loops with a higher migration energy when the number of 
atoms in the dislocation loop is greater than N;;,, and the corresponding value is 0.90 eV [29]. 
Jansson et al. used N= 90 as the transition threshold, and the corresponding migration 
energy was 0.2 eV [29]. We selected the same values for these two effective parameters in this 


study. For interstitial clusters with N; > 7, the attempt frequency is calculated as follows: 


Vinig Ni) = Gene (13) 


where the exponential part, e, of N, varies theoretically from 0.50 to 1.00, corresponding to the 
SIA cluster diffusion model [30] and kink movement model of dislocation loop edges [31], 
respectively. In this study, the value is taken as 0.8, as in reference [29]. 

Similar to the definition of stability for vacancy-type clusters mentioned previously, the 
dissociation energy of an SIA cluster is the sum of the binding energy of an SIA to its mother 
cluster plus and the migration energy. For the binding energy of the SIA to its mother cluster, we 


use the fitting function given by Domain et al. [28]: 


Ei, (N;) = 4 + 5.107249[(N,)3 — (N; + 1)3]. (14) 

The final assumption of this model is the use of a simple stochastic description to simulate 

SIA defect traps. The simulation does not explicitly introduce trap objects but defines the 

probability of an SIA-type defect transforming into a pinned defect. The transition probability is 

calculated using the relevant sink strength formula, as in reference [32]. The main trap type 
considered in the iron matrix is a C-vacancy trap; thus, its concentration is set to 50 appm. 


4. Results and discussion 


4.1 Evolution of Cu clusters during thermal aging 
Based on this model, we simulated the thermal aging process of copper clusters under the 
same conditions as the experiments performed by Chen et al. [33] and Vincent et al. [34] in a 


simulation box with dimensions of 72 x 72 x 72 nm°. The experimental alloys with Cu contents 


of 0.5 at.% and 0.6 at.% were annealed at 450 °C and 500 °C, respectively, and the precipitation 
of Cu clusters was characterized by atom probe tomography (APT). 

Because the concentration of vacancies in the simulation box is much higher than that in the 
real physical environment, the simulation time must be rescaled before comparison with the 
experimental data on the same time scale. Generally, the ratio of the Monte Carlo simulation time, 


tmc, to real physical time, treal, is defined as follows [13]: 


Che (15) 


v 
Creal 


Creal = tmc 


where Cy.q) is the experimental vacancy concentration, and Cyc is the simulation vacancy 
concentration. Considering that the interaction between solute atoms and vacancies will change 
the local vacancy formation energy during the precipitation of solute clusters, the effective 
concentration of vacancies in the simulation needs to be corrected [35]. In this study, we 


approximate Cyc using the mobile concentration [14], Cy opije, of vacancies as follows: 


v — _Njump 
Cmobile = AtAVv(pty’ (16) 


where AV v(pt) is the sum of the capture volumes of movable particle objects at each moment (i.e., 
the vacancies in this case), and Njymp is the number of vacancy jumps in the Fe matrix in time 
interval At; this value can be output directly from the simulation in MMonCa. It is worth 
mentioning that its ratio to the absolute vacancy concentration, fy, is positively related to the sink 
strength of Cu clusters to vacancies. 

For the real experimental vacancy concentration, the equilibrium estimation in a pure Fe 
environment can be obtained as follows: 

—hy 
Chea) = exp (op) (17) 
where hf is the enthalpy of formation of the vacancy in a pure Fe environment, which is taken to 
be 1.7 eV. 

Figure 3 shows a comparison of the thermal aging of the Fe-Cu alloy under different 
conditions, in which the OKMC simulation results are plotted as lines and correspond to the 
experimental results plotted as triangles. It is clear that our model is consistent with the 
experimental results and successfully shows the evolution of the density and size of the clusters. 

At 450 °C, as shown in Fig. 3(a) and (b), the average number density and radius of clusters 


versus time in the Fe-0.5 at.% Cu alloy increase gradually and then reach a steady state after 200 


h; this behavior is closely reproduced at the temperature investigated. The APT results show that 
the micro-hardness reaches maximum after 200 h of annealing, which may indicate the largest 
value of the cluster number density multiplied by the average radius [33]. 

At 500 °C, the number density increases rapidly at the initial stage of thermal aging and then 
gradually decreases, indicating obvious coalescence between the clusters, as shown in Fig. 3(c). 
The average radius of copper clusters was also measured during thermal aging at 500 °C for up to 
20 h, from which three different cluster growth processes were observed. In the first nucleation 
stage, the continuous precipitation of small-sized clusters is observed. Then, in the growth stage 


1/2 


of the clusters, the average growth rate of the clusters is close to the classical time law t “4^. Finally, 


the growth of the average cluster size is further reduced and approaches the classical time law 7', 
corresponding to the coarsening stage [36]. Simultaneously, the number density of clusters 
decreases in the later stage of cluster growth. For this result, the prediction of the OKMC model is 
in good agreement with the experimental data. 
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Fig. 3 (Color online) Average cluster radius and average cluster number density in the alloys 


after annealing at 450 °C and 500 °C 


Through comparison with the experimental results, we verified the validity of the model in 
this study. Moreover, to explain the evolution mechanism of copper clusters under different 
conditions, the parameters of the simulation process are further analyzed. As mentioned previously, 
in the simulated annealing process, we defined a recalculation for the simulation time; then, there 
is a fraction, fy, of the vacancy concentration between the mobile concentration and absolute 
concentration, which should be a variable of temperature and time. The changes in fy over the 
Monte Carlo time are shown in Fig. 4. It is clear that the fraction decreases with increasing 
simulation time, which indicates that the growing Cu clusters become a stronger sink for vacancies, 
and vacancies will then spend more time inside the clusters. It is also noticeable that there is a 
minimum value over the curing period at 500 °C, which may be related to the decrease in the 


average cluster number density according to the evolution of the copper clusters shown above. 
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Fig. 4 (Color online) Fraction of time spent by a vacancy in a pure iron environment at 450 °C 
(blue) and 500 °C (orange) 

Furthermore, the cluster density during thermal aging decreases because of the merging of 
large clusters and cluster dissociation through the emission of Cu-vacancy pairs, which can be 
observed by monitoring the simulation process. The one-by-one emission of Cu-vacancy pairs 
from mother Cu-vacancy clusters is important for the formation of large Cu clusters; when the Cu 
clusters maintain proper size and mobility, the merging of large clusters is possible. Considering 
the temperature effect, the Cu cluster grows much faster at 500 °C than at 450 °C. Because the 


cluster mobility is highly related to temperature, clusters have a better chance of merging together 


at higher temperatures. At the same time, it is assumed that the emission rate of Cu-vacancy 
couples reaches a balance with the absorption rate, after which the entire microstructure evolution 
process appears very slow; thus, there some of Cu atoms remain dispersed in the matrix. 
4.2 Evolution of Cu clusters under irradiation 

In this section, OKMC simulations were conducted for neutron irradiation under the 
representative conditions reported by Liu et al. [37] and ion irradiation experiments by Shu et al. 
[38]; the simulations were performed in a simulation box with dimensions of 72 x 72 x 72 nm. 

Chen33] obtained the size distribution of Cu clusters using histogram statistics, and the results 
showed a unimodal distribution with a peak at approximately 1.6 nm. In the neutron irradiation 
experiment performed by Liu[37], the peak size was approximately 1.45 nm. We carried out a 
neutron irradiation simulation in Fe-0.5 at.% Cu and compared the results with those from the 
annealing simulation, as shown in Fig. 5. It should be noted that the kernel density estimation 
(KDE) is used to obtain the cluster size distribution instead of a histogram because of the small 
sample size in the simulation box; correspondingly, its ordinate values are the probability density. 
This approach provides an estimation with excellent reliability, which can quickly stabilize with 
increasing samples by repeating the simulation. The results show that the size distribution obtained 
by the simulation of the thermal aging process has a peak value of 1.6 nm, whereas for the neutron 
irradiation simulation, the peak value is approximately 1.4 nm. This result is consistent with the 
experimental results, indicating that the cluster segregation produced by annealing at 450 °C for 


200 h is comparable to that under neutron irradiation at 290 °C with a flux of 5.4x10!8 n/cm?. 
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Fig. 5 (Color online) Size distributions of Cu clusters after irradiation at 290 °C and annealing at 
450 °C 

Furthermore, we compared the neutron irradiation and ion irradiation in Fe-0.8Cu using a 
damage rate of 3x10° dpa/s to simulate ion irradiation and a damage rate of 3x107 dpa/s to 
simulate neutron irradiation; the simulation results are shown in Fig. 6. 

The simulation results for neutron irradiation at lower dpa levels are shown as black lines on 
the left, corresponding to the experiments plotted as blue dots. The simulation results for ion 
irradiation at higher dpa levels are shown as red lines and correspond to the ion irradiation 
experiments plotted as purple triangles. Under irradiation, the number density of Cu-vacancy 
clusters formed in the alloy initially increases and then decreases, accompanied by an increase in 
the mean radius of the clusters, which is similar to the annealing simulation and corresponds to 
nucleation, growth, and coarsening stages. Meanwhile, the mean radius and number density ratios 
of precipitates in Fe-Cu alloys under ion irradiation are much smaller and larger, respectively, 


compared with those under neutron irradiation. 
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Fig. 6 (Color online) Average radius and number density of clusters in the alloy under neutron 
irradiation and ion irradiation at 290 °C 

For the irradiation simulations, it is found that under the same irradiation damage dose, the 
mean radius and number density of clusters under neutron and ion irradiation are significantly 
different; thus, we further analyze the vacancy clusters in the solute clusters. Figure 7 shows the 
2D distribution of copper and vacancies in solute clusters; it should be noted that this distribution 
is also obtained by the 2D KDE. Under neutron irradiation, from 0.0001dpa to 0.001dpa, the 
increase in the number of vacancies is small, as shown in Fig. 7(a) and (b). Thus, we can infer that 
the cluster maintains a high mobility, and the cluster merger mechanism occupies the later stage 
of evolution. The average radius of the clusters still increases rapidly, as shown in Fig. 6. Under 
ion irradiation, the content of vacancies in the cluster is high even at 0.03 dpa (Fig. 7(c)), which 
also leads to weakening of the cluster mobility. It can be inferred that the absorption of Cu-vacancy 


pairs should be dominant in the later stage of copper cluster precipitation such that the cluster size 
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reaches equilibrium, and most of the clusters contain less than 1000 Cu atoms, as shown in Fig. 


7(d). 
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Fig. 7 (Color online) Elemental distribution in copper—vacancy clusters at 290 °C under neutron 
irradiation at (a) 0.0001 dpa and (b) 0.001 dpa and ion irradiation at (c) 0.03 dpa and (d) 0.3 dpa 
Moreover, we compared the distributions of the binding energy, Egina Ncw Ny ), as shown in 

Fig. 8. It is clear that with an increase in N, in the Cu-vacancy cluster, the binding energy, 
Evina (Ncw Ny ), gradually increases. Under neutron irradiation, because of the slow damage rate, 
the vacancy accumulation rate is close to the vacancy dissociation rate in the cluster; therefore, the 
distribution in Fig. 7(b) is very similar to this binding energy distribution map. However, in the 
lower area of the map, a mutation range is observed. Therefore, when the irradiation rate is fast, 
i.e., corresponding to ion irradiation, the rapid accumulation of vacancies will lead to a rapid 
increase in the stability of the cluster. A similar tendency can be observed in the distribution shown 


in Fig. 7(c). 
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Fig. 8 (Color online) Binding energy of one vacancy to copper-vacancy clusters 

5 Conclusion 

In this study, we establish an OKMC model under the MMonca framework to simulate the 
solute precipitation process in Fe-Cu systems. The proposed model can be used to simulate the 
formation process of Cu-vacancy composite clusters under neutron/ion irradiation conditions, in 
addition to enabling prediction of the annealing process. Based on a modified description of the 
reaction network in the iron matrix, this model enables the simulation of sufficiently large volumes 
with long-term irradiation. The corresponding energies used for the parameterization are calculated 
from published DFT and AKMC results and provide a complementary simulation result to APT 
experiments. Our model includes all of the important mechanisms driving Cu precipitation in Fe. 

It is demonstrated that the merging of mobile Cu clusters containing thousands of Cu atoms 
is the main mechanism responsible for the growth and coarsening of Cu clusters during thermal 
aging. A more complex reaction network will occur in the Fe-Cu system under irradiation than 
under thermal aging, but there is still a mapping relationship in the cluster size distribution between 


both processes. In addition, the evolution behavior of the copper clusters differs under neutron and 


ion irradiation. The thermal aging results confirm that the merging of clusters is an important 
mechanism for the growth of copper clusters. However, because of the strong reaction between 
vacancy clusters and copper clusters, both clusters easily become nucleation sites for each other 
under ion irradiation, thus inhibiting the migration and coarsening of the copper clusters, ultimately 
leading to the smaller average size of copper clusters observed by APT experiments. 
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